---
title: "Van Egeren et al (2018) Fig. 4"
output: html_notebook
---


```{r echo=FALSE}
library(readr)
library(ggplot2)
library(reshape2)
library(RColorBrewer)
source("./utils.R")
```

## Fig. 4a-b

```{r include=FALSE}
filepath = "./"
dat = data.frame(A=numeric(0), B=numeric(0), C=numeric(0), D=numeric(0))
for (i in 2:6){
  for (j in c(0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99)){
    filename = paste(filepath, "data_var_r1/", j, "/", i, "/heritable/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel_prob = mean(temp$X2)
    filename = paste(filepath, "data_var_r1/", j, "/", i, "/heritable/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix_prob = mean(temp$X2)
    dat = rbind(dat, data.frame(A=i, B=tunnel_prob, C=fix_prob, D=j))
  }
}
colnames(dat) <- c("var", "tunnel", "fix", "r1")
dat["tunnel_norm"] = dat$tunnel/dat$fix

dat_simple = data.frame(B=numeric(0), C=numeric(0), D=numeric(0))
  for (j in c(0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 0.99)){
    filename = paste(filepath, "data_var_r1/", j, "/simple/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel_prob = mean(temp$X2)
    filename = paste(filepath, "data_var_r1/", j, "/simple/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix_prob = mean(temp$X2)
    dat_simple = rbind(dat_simple, data.frame(B=tunnel_prob, C=fix_prob, D=j))
  }
colnames(dat_simple) <- c("tunnel", "fix", "r1")
dat_simple["tunnel_norm"] = dat_simple$tunnel/dat_simple$fix
dat_simple["tunnel_SEM"] = sqrt((dat_simple$tunnel_norm) * (1-dat_simple$tunnel_norm) / (dat_simple$fix * 10000))
dat_simple["fix_SEM"] = sqrt((dat_simple$fix) * (1-dat_simple$fix) / (10000))

dat["tunnel_SEM"] = sqrt((dat$tunnel_norm) * (1-dat$tunnel_norm) / (dat$fix * 10000))
dat["fix_SEM"] = sqrt((dat$fix) * (1-dat$fix) / (10000))
```

```{r}
# Fig. 4b
ggplot(data=dat) + geom_line(data=dat_simple, aes(x=r1,y=tunnel_norm), size=1, color="#707070") + geom_errorbar(data=dat_simple, aes(x=r1,ymin=tunnel_norm-1.96*tunnel_SEM, ymax=tunnel_norm+1.96*tunnel_SEM),width=0, color="#707070") + geom_line(aes(x=r1,y=tunnel_norm,group=var,color=-var)) + geom_errorbar(aes(x=r1,ymin=pmax(tunnel_norm-1.96*tunnel_SEM,0), ymax=tunnel_norm+1.96*tunnel_SEM,color=-var),width=0) + theme_bw()+ scale_color_gradient2(low="#cccc00", mid="#00cc00", high="#0066ff", midpoint=-4) + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank()) + xlim(0.3, 1.03)

# Fig. 4a
ggplot(data=dat) + geom_line(data=dat_simple, aes(x=r1,y=fix), size=1, color="#707070") + geom_errorbar(data=dat_simple, aes(x=r1,ymin=fix-1.96*fix_SEM, ymax=fix+1.96*fix_SEM),width=0, color="#707070") + geom_line(aes(x=r1,y=fix,group=var,color=-var)) + geom_errorbar(aes(x=r1,ymin=pmax(fix-1.96*fix_SEM,0), ymax=fix+1.96*fix_SEM,color=-var),width=0) + theme_bw()+ scale_color_gradient2(low="#cccc00", mid="#00cc00", high="#0066ff", midpoint=-4) + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank()) 
```

## Fig. 4c

```{r}
newmut_dat = data.frame(fit=numeric(0),fix_density=numeric(0),total_density=numeric(0),lifetimes=character(0))
for (i in c("1", "5", "10", "20", "50", "100", "heritable")){
  filename = paste(filepath, "new_mut_processed/processed_", i, ".csv", sep="")
  temp = read_csv(filename, col_names = FALSE)
  colnames(temp) = c("fit","fix")
  #print(sum(temp$fix))
  #print(length(temp$fix))
  fix_density = density(temp[temp$fix==1,]$fit, bw=0.03, from=0.93, to=1.07, n=1000)
  total_density = density(temp$fit, bw=0.03, from=0.93, to=1.07, n=1000)
  newmut_dat = rbind(newmut_dat, data.frame(fit=fix_density$x, fix_density= (sum(temp$fix)/length(temp$fix)) * fix_density$y, total_density=total_density$y, lifetimes=i))
}
newmut_dat["cond_density"] = newmut_dat$fix_density/newmut_dat$total_density
```

```{r}
# diagnostic plot showing the smoothed density of the fitness of new intermediate stage mutants
#ggplot(data=newmut_dat[newmut_dat$lifetimes != "heritable",]) + geom_line(aes(x=fit, y=fix_density, color=sapply(lifetimes,to_colorscale), group=sapply(lifetimes,to_colorscale))) + theme_bw() + geom_line(data=newmut_dat[newmut_dat$lifetimes == "heritable",], aes(x=fit,y=fix_density), color="#0000cc", size=1) + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank()) + scale_color_gradient(low="#ccccff", high="#0000ff")

# Fig. 4c
ggplot(data=newmut_dat[newmut_dat$lifetimes != "heritable",]) + geom_line(aes(x=fit, y=cond_density, color=sapply(lifetimes,to_colorscale), group=sapply(lifetimes,to_colorscale))) + theme_bw() + geom_line(data=newmut_dat[newmut_dat$lifetimes == "heritable",], aes(x=fit,y=cond_density), color="#4400cc", size=1) + theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank()) + scale_color_gradient2(low="#cccc00", mid="#ff3399", high="#9900ff", midpoint=4)
```

## Fig. 4d-e

```{r include=FALSE}
# THIS WILL LIKELY TAKE A WHILE
variances = extract_vars(filepath)
her_variances = extract_vars_heritable(filepath)
```

```{r include=FALSE}
fix_dat = data.frame(A=numeric(0), B=numeric(0), C=numeric(0), D=numeric(0))
num_runs = 10000
for (i in 2:6){
  for (j in c("1","5","10","50","100")){
    filename = paste(filepath, "var_tau_sweep/", j, "/", i, "/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix = sum(temp$X2)/num_runs
    filename = paste(filepath, "var_tau_sweep/", j, "/", i, "/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel = sum(temp$X2)/num_runs
    tunnel = tunnel/fix
    fix_dat = rbind(fix_dat, data.frame(A=i,B=j,C=fix,D=tunnel))
  }
}
for (i in 2:6){
  for (j in c("heritable")){
    filename = paste(filepath, "var_tau_sweep/", j, "/", i, "/iftype2.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    fix = sum(temp$X2)/num_runs
    filename = paste(filepath, "var_tau_sweep/", j, "/", i, "/type_1_tunnel.oevo", sep="")
    temp = read_csv(filename, col_names=FALSE)
    tunnel = sum(temp$X2)/num_runs
    tunnel = tunnel/fix
    fix_dat = rbind(fix_dat, data.frame(A=i,B=j,C=fix,D=tunnel))
  }
}
filename = paste(filepath, "var_tau_sweep/simple/iftype2.oevo", sep="")
temp = read_csv(filename, col_names=FALSE)
fix = sum(temp$X2)/num_runs
filename = paste(filepath, "var_tau_sweep/simple/type_1_tunnel.oevo", sep="")
temp = read_csv(filename, col_names=FALSE)
tunnel = sum(temp$X2)/num_runs
tunnel = tunnel/fix
fix_SEM_hold = sqrt(fix * (1-fix)/num_runs)
tunnel_SEM_hold = sqrt(tunnel * (1-tunnel)/(fix * num_runs))
fix_hold = fix
tunnel_hold = tunnel
colnames(fix_dat) = c("input_var", "lifetimes", "fix", "tunnel")
fix_dat["log_var"] = c(variances$log_var, her_variances$log_var)
fix_dat["log_upper"] = c(variances$log_upper, her_variances$log_upper)
fix_dat["log_lower"] = c(variances$log_lower, her_variances$log_lower)
fix_dat["fix_SEM"] = sqrt(fix_dat$fix * (1-fix_dat$fix)/num_runs)
fix_dat["tunnel_SEM"] = sqrt(fix_dat$tunnel * (1-fix_dat$tunnel)/(fix_dat$fix * num_runs))
```

```{r echo=FALSE}
# Fig. 4d
ggplot(data=fix_dat[fix_dat$lifetimes != "heritable",]) + geom_hline(aes(yintercept=fix_hold), color="#707070") + geom_line(aes(x=log_var,y=fix,group=lifetimes,color=sapply(lifetimes, to_colorscale)))+ geom_errorbarh(aes(x=log_var,xmin=log_lower,xmax=log_upper,y=fix, color=sapply(lifetimes, to_colorscale))) + theme_bw() + geom_errorbar(aes(x=log_var, ymin=fix-1.96*fix_SEM, ymax=fix+1.96*fix_SEM, color=sapply(lifetimes, to_colorscale)))+ theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank())  + scale_color_gradient2(low="#cccc00", mid="#ff3399", high="#9900ff", midpoint=4) + geom_line(data=fix_dat[fix_dat$lifetimes == "heritable",], aes(x=log_var,y=fix,group=lifetimes), color="#4400cc", size=1) + geom_errorbar(data=fix_dat[fix_dat$lifetimes == "heritable",], aes(x=log_var, ymin=fix-1.96*fix_SEM, ymax=fix+1.96*fix_SEM), color="#4400cc", width=0)+ geom_errorbarh(data=fix_dat[fix_dat$lifetimes == "heritable",], aes(x=log_var,xmin=log_lower,xmax=log_upper,y=fix),color="#4400cc", height=0)

# Fig. 4e
ggplot(data=fix_dat[fix_dat$lifetimes != "heritable",])  + geom_hline(aes(yintercept=tunnel_hold), color="#707070") + geom_line(aes(x=log_var,y=tunnel,group=lifetimes,color=sapply(lifetimes, to_colorscale)))+ geom_errorbarh(aes(x=log_var,xmin=log_lower,xmax=log_upper,y=tunnel, color=sapply(lifetimes, to_colorscale)), height=0) + theme_bw() + geom_errorbar(aes(x=log_var, ymin=tunnel-1.96*tunnel_SEM, ymax=tunnel+1.96*tunnel_SEM, color=sapply(lifetimes, to_colorscale)), width=0)+ theme(legend.position = "none", axis.title.x=element_blank(), axis.title.y=element_blank(), text = element_text(size=20), panel.grid.major = element_blank(), panel.grid.minor = element_blank(), axis.line = element_line(colour = "black"), panel.border = element_blank()) + scale_color_gradient2(low="#cccc00", mid="#ff3399", high="#9900ff", midpoint=4) + geom_line(data=fix_dat[fix_dat$lifetimes == "heritable",], aes(x=log_var,y=tunnel,group=lifetimes), color="#4400cc", size=1) + geom_errorbar(data=fix_dat[fix_dat$lifetimes == "heritable",], aes(x=log_var, ymin=tunnel-1.96*tunnel_SEM, ymax=tunnel+1.96*tunnel_SEM), color="#4400cc", width=0)+ geom_errorbarh(data=fix_dat[fix_dat$lifetimes == "heritable",], aes(x=log_var,xmin=log_lower,xmax=log_upper,y=tunnel),color="#4400cc", height=0)
```